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An initially homogeneous freely evolving fiuid of inelastic hard spheres develops inhomo- 
geneities in the flow fleld u(r, t) (vortices) and in the density field n(r, t) (clusters), driven 
by unstable fiuctuations, Sa = {Sn,Su.}. Their spatial correlations, ((5a(r, t)<5a(r', t)}, as 
measured in molecular dynamics simulations, exhibit long range correlations; the mean 
vortex diameter grows as ^(t) oc Vint; there occur transitions to macroscopic shearing 
states, etc. 

The Cahn-Hilliard theory of spinodal decomposition offers a qualitative understanding 
and quantitative estimates of the observed phenomena. When intrinsic length scales are 
of the order of the system size, effects of physical boundaries and periodic boundaries 
(finite size effects in simulations) are important. 

1. Introduction 

Standard fluids, when out of equilibrium, wifl rapidly decay to local equilibrium 
within a few mean free times. The subsequent decay of spatial inhomogeneities to- 
wards global equilibrium is controlled by the slow hydrodynamic time evolution. A 
prototypical model for this is a system of N smooth elastic hard spheres. Compare 
this system to a system of smooth inelastic hard spheres, defined such that in the 
center of mass frame a fraction of order e — \ ~ a (where a is called the restitution 
coefficient; < a < 1) of the kinetic energy is lost in each collision (for precise 
definitions we refer to Ref. 1). Linear momentum is conserved during collisions. 
Consequently the local momentum density, g(r, i) = n{r,t)u{r,t), or the flow ve- 
locity u(r, t) is a slowly varying quantity, and the system can still be considered a 
fluid. The lack of energy conservation makes this fluid, whether driven by gravitv 
or shear stresses, or freely evolving, behave very differently from a standard fluid.a 
The system of inelastic hard spheres represents an idealized model for rapid 
granular flows, where the dynamics of individual (macroscopic) particles is described 
by binary collisions, separated by free propagation over a typical mean free path. 
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To study this problem by computer simulations, we consider a 2-dimcnsional 
system of N inelastic hard disks of diameter a and unit mass, contained in a volume 
V = L^ with neriodic boundary conditions, and use an event-driven molecular 
dynamics code.u At the initial time the system starts off in a spatially homogeneous 
equilibrium state of elastic hard disks (e = 0). From there on the kinetic energy is 
dissipated by collisions, which leads to cooling. 

We monitor the time evolution of the system, and measure the velocity dis- 
tribution of a particle, the time dependence of the kinetic energy, and the total 
number of collisions. These are the standard quantities that can be obtained from 
the Enskog-Boltzmann equation. 

The novel feature of the present study is an analysis of the spatial correlations 
in mass and energy densities, and in flow velocities at different points in the system, 
which exhibit long range spatial correlations. To understand such dynamic corre- 
lations one has to go beyond the standard Enskog-Boltzmann equation, which is 
based on the molecular chaos assumption, and use the Ring kinetic theory, which ac- 
counts for sequences of correlated binary collisions. Ring kinetic theory and closely 
related mode coupling theories have been succesful in explaining the phenomena of 
long time tails, and of long range spatial correlations in nonequilibrium stationary 
states,QEl and in other systems that violate the conditions of detailed balance.l3 This 
violation also occurs in systems of inelastic hard spheres. 

In fact, we have applied the Ring theory to such systems, and succesfuUy ex- 
plained the behavior of the nonequilibrium pair distribution at late times and large 
distances. Preliminary results of this investigation have been reported at several 
conferences and workshops. 

The plan of this paper is as follows: section 2 describes the homogeneous cooling 
state. Spatial fluctuations (section 3) drive the system away from this state through 
slow hydrodynamic modes (section 4). Some modes are unstable and lead to phase 
separation and clustering (section 5). Vorticity diffusion controls the long time 
behavior (section 6) and leads to scaling laws. Section 7 deals with the effects of 
boundaries on the structure of asymptotic states. 

2. Homogeneous Cooling State 

What are the important observations that can be deduced from molecular dynamics 
simulations, performed on a system of N inelastic hard disks, that are prepared in 
a spatially homogeneous equilibrium state of an elastic hard disk system? 

The system stays for many collision times in the so-called homogeneous cool- 
ing state, described by a single particle distribution function which is essentially a 
Maxwellian with a time dependent total kinetic energy A^£'(i).aa'El As long as the 
system is in this homogeneous cooling state, the energy per particle equals the tem- 
perature, T{t) = E{t). This temperature decays as t~^ for long times (see Fig. 1), in 
quantitative agreement with the predictions from kinetic theory, and has the form 

T{t) EE lv^,it)=TiO)/[l + t/t,r 



Long Range Correlations in Idealized Granular Flows 3 



= T(0)exp[-27oT]. 



(1) 



Here 70 = |(1 — a^) = ie(l — ^e) measures the degree of inelasticity, and vo{t) is 
the thermal velocity. The characteristic time of homogeneous cooling, te = io/70; 
is proportional to the mean free time for elastic hard disks in the initial state, 
to = l/aJo(0) = l/lV^nnxcrvoiO)], where x is the radial distribution function for 
elastic hard disks at contact. With this definition, the mean free path at the initial 
time equals that of the elastic hard disks and is given by ^0 — fo(0)io- The average 
number of collisions r, suffered by a single particle in a time t, follows by integrating 
dr = LOQ{t)dt, where wo(i) is the collision frequency of inelastic hard disks, which is 
proportional to y/T{t). The result is 

7oT = ln[l+t/g, 

where r is related to the total number of collisions C as C = \tN ■ 



(2) 
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Fig. 1. Loft: Plot of the logarithm of the total energy as a function of the number of collisions 
per particle t in a simulation with A'^ = 20000 particles, system size of L = 253cr, and dissipation 
parameter e = 1/40. The straight (dashed) line with slope —270 ~ —0.0247 (Eq. (1)) corresponds 
to the behavior of the homogeneous cooling state. Deviations of E{t) from the straight line 
around t = 200 — 250 are correlated with the appearance of inhomogeneities in the density field 
(see behavior of Gnn for r = 170 — 332 in Fig. 2). Right: Number of collisions per particle 
T versus time t for the same system. The straight dashed line correponds to the homogeneous 
cooling prediction, Eq. (2). 



In the subsequent time evolution of the inelastic hard disk system, the homoge- 
neous cooling state plays a role very similar to local equilibrium in ordinary fluids, 
and provides a conceptual basis for a hydrodynamic description. However, this ho- 
mogeneous cooling state does not at all behave as an equilibrium state of TV elastic 
hard disks with a time dependent temperature. In the latter state spatial correla- 
tions are caused by hard core excluded volume effects and extend only over a few 
disk diameters a. In the present system of inelastic hard spheres, however, we have 
observed: (i) long range spatial correlations extending over rnore than a decade of 
disk diameters; {ii) slowly growing spatial inhomogeneitiesav^ in the local flow field 
u(r, i) and density field n{r,t), which can be described by coupling of hydrody- 
namic modes, at least for times that are not too long; (Hi) a collision number t that 
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increases^ more rapidly than (I/70) ln(l + t/te), and a cooling law that decays more 
sZow/ipllj than the prediction of Eq. ffl), which is iUustratcd in Fig. 1. 

3. Spatial Correlations 

To obtain a better understanding of the dynamics of a freely evolving system of these 
inelastic hard disks, we follow the local mass density n(r,t), the local momentum 
density g(r, t) — n(r, i)u(r, t) and the local energy per particle e(r, t) in a single 
realization of a system of A^ = 20000 particles by means of computer simulations, 
where the inelasticity parameter is e = 1/40 and the volume fraction (p — ■^nna'^ ~ 
0.245, corresponding to a system size L ~ 253(7. 

Visual observation of snapshots of n{r,t)- and g(r, t)-fields show that at early 
times (r = 82), the momentum density does not exhibit much visible structure, and 
the mass density is totally homogeneous. At r = 120— 170 the mass density is still 
homogeneous, but the flow field starts to show vortices after coarse graining over 
cells of 5(7 X 5(7. At the latest time (r = 332) the flow field shows a pronounced 
structure of vortices, also in the fine grained u(r, i)-field, with a typical diameter ^ 
of the order of j of the system size L, and the mass density starts to show barely 
visible inhomogeneities in the fine grained n(r, i)-field, but clearly visible in the 
coarse grained n(r, i)-fields. 

To quantify these observations we have measured the spatial correlation func- 
tions of the microscopic densities, i.e. 

7i^Gab{r,t) = ^ f dR{5a{R + r,t)6b{R,t)). (3) 

Here the labels a, b take the values {n, e,a, ||,_L}, and refer respectively to the 
microscopic mass density, momentum density and local energy per particle, 

Sfi{r,t) = ^(5(r,(t)-r)-n 

i 

niLa{r,t) = ^Via{t)S{ri{t) - r) 

i 

nSe{r,t) = ^ i (^K*) - ^o'W) '^(r.W - r). (4) 



The carets denote microscopic quantities. Moreover, vo{t) — y/2T{t) is the thermal 
velocity and v^, r^ are the velocity and position of the i-th particle. Greek labels 
a, (3 — {x, y} refer to Cartesian components. The second rank tensor field Gapi'r, t) 
has transverse (_L) and longitudinal (||) components 

G||(r, t) = TarpGa/sir^t) 

G±{r,t) ^ f_Laf±/3G'Q/3(r,t), (5) 

where f , f^ are unit vectors, parallel and perpendicular to the relative position r. 
As long as the system remains spatially homogeneous with a vanishing flow velocity, 
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the local energy and temperature fluctuations are the same. However, as soon as 
the average flow velocity u(r, t) is nonvanishing, the granular temperature is the 
average energy per particle in the local rest frame. 

The measurement of correlation functions in a single realization is shown in Fig. 
2. The results are remarkable! Even at the earliest time (r = 82), long before 
any vortices are visible, long range spatial correlations have developed between the 
local flow velocities at different points in the fluid with a typical correlation length 
^ ~ 30(7 (f is the location of the minimum in G±{r,t)), which gradually increases 
to about 60cr at r = 332. The negative minimum in G±{r,t) clearly signals the 
presence of vortices with a diameter ^. It is also interesting to compare the noise 
level and magnitude (vertical axes in Fig. 2) of all fluctuations, relative to the 
homogeneous cooling state, and correlate the signals above noise level and their 
growth rates with the visual observations of the n(r, t)- and g(r, i)-snapshots. 

The dynamics of the phenomena observed seems to be controlled for a long time 
by the transverse flow field or vorticity modes. We conjecture that there exists an 
extended time regime in which these modes describe the behavior of the spatial 
correlation functions, at least of those correlation functions that involve the flow 
fields. We shall test this conjecture by analyzing the simulation results. 

Moreover, we note that in normal fluids with elastic collisions, times t > 5 
are already considered long in kinetic theory. The long time tail in the velocity 
autocorrelation function of elastic hard spheres is typically observed in the ranges 
10 ^ ''' ^ 60. The range of the spatial correlation functions G±{r,t) and G'||(r, i), 
observed here, extends far beyond any static correlation length (a few cr's) and far 
beyond the mean free path (Zq ~ 0.8cr at this density). 

4. Hydrodynamic Dispersion Relations 

The previous observations on large spatial and temporal scales suggest that the long 
wavelength components of the microscopic fluctuations are described by linearized 
(because we are dealing with fluctuations) hydrodynamics. This should hold as long 
as fluctuations have not grown to macroscopic size. Beyond this time the evolution 
is controlled by the full nonlinear hydrodynamic equations. 

The fluctuations discussed in the previous sections are taken relative to the ho- 
mogeneous cooling state. Consequently, we linearize the hydrodynamic equations 
around this state,aEJ and consider the Fourier components of the hydrodynamic 
fields 5n(k, t), Ua{'k,t)/vo{t) and 5e(k, t)/w§(i) where T(i) = ^Wo(0 is the tempera- 
ture in the homogeneous cooling state. The resulting set of linear equations contains 
still time dependent coefficients, such as the local pressure, proportional to T(i), 
and transport coefficients, which are assumed to be independent of the dissipation 
parameter e, and proportional to the average collision frequency ujQ{t) oc ^jT{t). 
Hence, the kinematic viscosity of inelastic hard disks is v{t) = z/exp(— 7or), where 
v — rj/mn is the corresponding viscosity of elastic hard disks at the initial time. 

By changing to the time variable r, introduced in Eq. (y), the linearized hydro- 
dynamic equations transform into a set of coupled linear equations with constant 
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Fig. 2. Spatial correlations of relative fluctuations around the homogeneous cooling state 
Gabi^^t)/^^ in mass densities (nn), flow velocities (||, -L) and energy densities (ee) (p = 0,1,2 
respectively), plotted versus r/cr at a fixed number of collisions r = 2C'/N per particle, measured 
in a single run on a system with TV = 20000 particles in a volume V = L^ with L = 253cr and 
dissipation parameter e = 1/40. Notice the vertical scales, showing magnitudes and growth rate 
of the fluctuations: at r = 332, Gii and G± have increased (20x), Gnn (lOx) and Gee (Ix)- 
The measurements of GQi,(r,t), shown here, are coarse grained by collecting data points in bins of 
width Ar = 5cr. For bin sizes of Ar = 0.2(7 there is already a 'signal over noise level' in Gii and 
Gx at T = 82, whereas the signals in Gnn and Gee just start to rise above noise level at short 
distances r < 30cr at r ~ 220. 
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coefficients. The resulting eigenmodes have the form ■(/'/j(k, t) = ipf^{'k)exp[z^{k)T], 
where /i — {J-,H,(t — ±} labels the shear mode (_L), the heat mode (H) and the 
sound modes (ct = ±). Note that the modes only show exponential behavior if time 
is measured by counting collisions. In real time the modes decay algebraically. The 
shear mode ■(/'-L(k, r) is u±(k,t)/vQ(t), and the heat and sound modes are linear 
combinations of the relative fluctuations Jn(k, i), u\\{'k,t)/vo{t) and (5e(k, i)/wg(t). 
Here the subscripts longitudinal (||) and transverse (_L) refer to the direction of k. 
The dispersion relations for the eigenvalues z^ of heat and shear mode are shown 
in Fig. 3 as a function of the rescaled wavenumber q = ka/ ^/jq. These functions 
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Fig. 3. Dispersion relations for the decay rates 2,j of hydrodynamic modes {fi = {±, H, it}) in units 
70, plotted versus the reduced wavenumber q = ka/^/^ to make z±/^o independent of dissipation 
parameter 70 = 9^(1 ~ o^)' Numerical values are calculated from the Enskog-Boltzmann equation 
for elastic hard disks at volume fraction <f> = jnna^ ~ 0.245. Wavenumber q'^ and g^ are the 
thresholds of stability for heat and shear modes. 



intersect the q-axis respectively at q'^ and g^, where 5* = fc^cr/-\/7o- At volume 
fraction (f> ~ 0.245 one has q'^ ~ 0.50 and q]_ ~ 1.17. The numerical values of 
the transport coefficients are calculated from the Enskog theory for elastic hard 
diskgl3 at a volume fraction cj) ~ 0.245 and initial temperature T(0) — 1, used in 
the present simulations. We only quote the result for the shear mode or transverse 
velocity field ,aElO which decays as 



uj_{'k,t) ^ vo{t)c^p[z±{k)T] 

= «o(0)(l+tAe)-''''*= 



(6) 



where the dispersion relation for the relaxation rate is: 

z±{k) = 70(1 - vk'^te) = 7o(l - i>q'^), 



(7) 



with i> = vto/a'^. If the viscosity of the inelastic hard disks depend on the dissipation 
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parameter e only through ^jT{t) — as we have assumed here — , then zx(fc)/7o is a 
universal function of g = ka / ^FJq^ independent of the dissipation parameter 70. 

5. Instabilities and Clustering 

The spectrum of eigenvalues in Fig. 3 shows that the relative fluctuations with 
respect to the homogeneous cooling state contain several unstable modes with an 
exponential growth rate, exp[z^(fc)r], if time is measured in collision numbers r. 
These unstable modes, ■i/'/j(k, r), are the long wavelength shear mode (/i =_L) with 
k < k\, and the long wavelength heat mode ipnik, t) with k < k*^ < k\. The sound 
modes are stable for all wavelengths. The dynamics of these slow modes determines 
the time evolution of the long wavelength fluctuations, away from the homogeneous 
cooling state, and determines a mechanism for initial pattern selection. In the initial 
stages spontaneous fluctuations with respect to the homogeneous cooling state occur 
on all possible wavelengths. Those with k smaller than k*jj or k\_ will grow at an 
exponential rate exp[2;^(fc)T]. 

The fastest growing mode is the one with the smallest wavenumber, as can be 
seen from the dispersion relations in Fig. 3. In finite systems with periodic boundary 
conditions, as used in the present computer simulations, the smallest wavenumber 
allowed is k„i — I-k jL. As long as fc„i < /c|^, unstable shear modes as well as heat 
modes can be excited, and we can estimate the onset times t^^ for the shear and 
heat mode instabilities by the criterion 

z^(fc„)r^~l (^=±,H), (8) 

which implies that the amplitude of relative fluctuation has increased by a factor 
e above its normal level. As long as the fluctuations are small and, hence, can be 
described by linearized hydrodynamic equations, the vorticity mode does not couple 
to density and energy fluctuations, and the onset for the appearance of vortex 
structure is the time tj_ , as defined in Eq. (H) . The fluctuations in mass and energy 
densities, however, do couple to the unstable heat mode, and the onset time for 
observable inhomogeneities in the mass density (clustering) and energy field will be 
of order th- 

For the system used in the simulations of section 3 with L — 253(T and e = 0.025 
one has qm — 0.22. The corresponding rate constants Zf^{km) can be read off from 
Fig. 3 and yield estimates of the onset times for the appearance of vorticity {t± ~ 83) 
and for clustering in the density {th — 225). Inspection of the configurations n(r, t) 
and g(r, t), as well as the correlations in Fig. 2, shows that the structures in Gnnir, t) 
and Gee{r,t) are barely noticeable at r = 170, and well developed at t = 332, 
consistent with the theoretical estimates for t±_ and th- 

Suppose we decrease the system size or decrease e, such that qm = 2Tra/[Ly/jQ] 
increases and moves into the interval (g|^,g^). Then only shear modes can be 
excited and the onset time t±_ for the appearance of vortex structures will increase 
like l/z±{km) with increasing q^- The vorticity will grow, but density and energy 
fields will remain spatially homogeneous, as long as linearized hydrodynamics is 
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adequate. An example of such a system is realized at iV = 5000, L = 1260" and 
e — 0.007, where g™ ~ 0.84, which is located in the middle of the interval {q'^, q'^). 
The simulations show that the total energy follows the homogeneous cooling law, 
E{t) ~ i?(0) exp[— 27ot], up to t ~ 1100. Of course, in all cases where km < fc^, 
eventually the fluctuations will growto macroscopic size, and after a very long time 
nonlinear hydrodynamics takes over ,□ and leads again to inhomogeneities in density- 
and energy fields. 

If one still further increases g,„ , it will pass the stability threshold g^ , and none 
of the unstable modes is excitable. All hydrodynamic modes become stable, and the 
system will remain in the homogeneous cooling state for all times. McNamara and 
YoungtJ call this state the kinetic state. An example is the system {N ~ 5000, L = 
126a, e = 0.001}, where qm — 2.2 > q*^. Even after t — 2044 no inhomogeneities 
in density or flow fields have developed, and the homogeneous cooling law is still 
obeyed, as can be inferred from the measured correlation functions G'at(r, i) and 
the measured kinetic energy E{t). 

6. Vorticity Diffusion and Scaling 

A rough estimate of the behavior of the pair correlation function can be obtained 
from the Cahn-Hilliard theorytJ for the dynamics of phase separation. In this theory 
one calculates the structure factors, defined as the Fourier transforms S'ah(k, f) = 
V~^{5a{'k, t)Sb{'k,t)) of the correlation functions Gab(r, t), from the long wavelength 
behavior of the unstable hydrodynamic (macroscopic) modes (5a (k, t). Very recently 
Deltour and Barrat used this theory to calculate the structure factor S'„„ for a fluid 
of inelastic hard disks. Ej 

Here a similar estimate for the structure factor S±{k,t) = V^^{\u±(]s.,t)\'^) 
of the vorticity field yields S±{k,t) ex VQ{t) exp[2z±{k)T] ex exp[~2i'tok^T]. On 
the basis of the conjecture, formulated in section 3, the contributions of heat and 
sound modes can be neglected. The correlation functions G±{r,t) and G\i{r,t) 
can then be obtained by taking inverse Fourier transforms of (k • r)'^S±{k,t) and 
[1 — (k • r)'^]S±{k,t) respectively, where a = a/|a|. They have the generic scaling 
form T^^ F^(r / y/vtor) with /i = {_L, ||}. The explicit form of F^{x) is not used in 
the arguments presented below. 

We recall that the average diameter ^(t) of a vortex can be identified with the 
location of the minimum of G±{r, t). It then follows from the scaling form and Eq. 
(0) that 



^{t) ~ ao^utaT ex ^j/te ln(t/te), (9) 

where the constant of proportionality ao is independent of the volume fraction (j) 
and of the dissipation parameter e. For large systems the value of ao, measured in 
the simulations, is typically oq ~ 4.0 — 5.3 or 3.8 — 4.0 in the r-intervals (10 — 80) or 
(82 — 216) respectively for systems with parameters {N = 5 x 10*, e = 0.1, (/) = 0.4} 
or {A^ = 2 X 10"*, e = 0.025, (/) = 0.245}. For small systems (typically A^ = 5000, 
e < 0.1), ^(i) saturates to a constant value and diffusive growth is absent due to 



10 Long Range Correlations in Idealized Granular Flows 

interference effects caused by the periodic boundaries. 

The diffusive growth of the average vortex diameter offers strong support for 
the vahdity of our conjecture that the long time dynamics on large spatial scales is 
mainly governed by vorticity diffusion. More refined Ring kinetic theory or mode 
coupling theories provide detailed informationEj about the explicit analytic forms 
of Gfi{r,t). It is found that data collapse for G^{r,t) with /i = {_L, ||} only hold 
for large distances, where r > ^i/taT, and large times, where t > t^ — (ln2)/7o. 
Here Te = T(ie) is the characteristic time scale for homogeneous cooling, as defined 
through Eq. (|[). 

We conclude that our theoretical analysis has provided a consistent picture of 
the long range dynamic correlations that have built up through the dissipative 
dynamics of inelastic hard sphere systems. This dynamics violates the conditions 
of detailed balance. 

The theoretical analysis also predicts the existence of algebraic long range cor- 
relations oc 1 /?''*, where d is the dimensionality of the hard sphere system. Similar 
algebraic tails have been observed in lattice gas automata with collision rules that 
violate the conditions of detailed balanced The complete analysis and theoretical 
predictions, which have no adjustable parameters, as well as extensiva comparisons 
with the results of computer simulations will be reported elsewhere.tJ 

A comparison with the simulation results for a single realization is shown in 
Fig. 4. 
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Fig. 4. Comparison of theory (solid line) and simulations (circles) of Gi\{r,t) and G±{r,t) as 
functions of r/a at r = 2C/N ~ 40 collisions per particle in a system with N = 5000, L = 126(T, 
dissipation parameter e = 0.1, homogeneous cooling length le ~ 7cr, with unstable heat and 
shear modes [qm = 27Ta/[L^/y)] ~ 0.23 < q^), and estimated transition to the shearing state at 
rtr = L'^/lGiuto] ~ 325. 



There is a quantitative agreement between theory and simulations for all r- 
values, provided r is sufficiently small (r < 60 in Fig. 4) to avoid effects from 
nonlinear hydrodynamics, and provided system size L and/or dissipation parameter 
e are sufficiently large to avoid interference effects from periodic boundaries. 
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7. Transition to Shearing States/Boundary Effects 

What are the asymptotic states observed in computer simulations? In large but 
finite systems with L and e large enough (that is qm < q^), so that all unstable 
modes can be excited, the mean vortex diameter can not grow indefinitely. It will 
typically stop growing at a transition time Tt^, where ^(rtr) — ^L, and the periodic 
boundaries force a transition to an inhomogeneous shearing state. The dynamics of 
this transition, caused by fusion of vortices of equal sign into shearing layers, is well 
illustrated in the snapshot of the fiow field in Fig. 10b of Rcf. 10. In this case, the 
transition time rtr can be estimated from the growth law ^(t) ~ A^i^toT, and yields 
Ttr — i^/64z^toj which is proportional to the volume of the system, and independent 
of the dissipation parameter e. For the system with {N = 5000, L — 126, e = 0.1} 
the theoretical estimate gives rtr — 325, whereas the simulations yield rtr — 400±25. 
In the system of Rcf. 2 with {iV = 4 x 10"*, L ~ 793, e == 0.4} the above estimate 
gives Ttr — 525, whereas the value reported for the simulation is rtr = 20 /N ~ 600. 

However, for systems with {N = 5000, L ~ 126, e < 0.01}, and for the small 
systems of Ref. 10 with {N ^ 1024:, L == 57, e ~ 0.18 - 0.45} diffusive growth of 
the mean vortex diameter is absent due to interference effects resulting from the 
periodic boundaries, as the stability wavelength for the heat mode Ih = ^ir/k^ is 
of order L (see discussion below). As diffusive growth is absent, the transition time 
to the shearing state, Ttr = 2C/N ~ 1400 — 1600, observed in the simulations of 
Fig. 5, and of Ref. 10, has no connection to the theoretical estimate Ttr = L^ /GAiyto, 
valid for large systems. 




Fig. 5. Flow field, coarse grained over cells of 3cr x 3a, showing a transition at Ttr = 2C/N ~ 1600 
collisions per particle to a shearing state with a homogeneous density distribution in a small system 
with N = 5000, L = 126o", dissipation parameter e = 0.01, stability wavelength Ih — 126o", and 
homogeneous cooling length Ze — 160(t (see section 7). 



In the inelastic hard sphere system under discussion, there exists a number of 
intrinsic dynamical length scales: the mean free path Iq, the homogeneous cooling 
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length le = wo(0)ie = lo/lo', the mean vortex size ^(i); the stabihty wavelength, 
Ifj, = 27r/fc*, for the stability of the heat {fi = H) and shear mode (/i =-L). 

If any of these length scales become of the order of the system size L, the effects 
of the boundaries dominate the behavior of the system. If one is dealing with finite 
geometries, the effects of physical boundaries become important, and one needs to 
model particle-wall interactions. 

In case one is dealing with computer simulations, in which periodic boundary 
conditions are used, the intrinsic length scales should be less than ^L, say, in order 
to minimize finite size effects. If, however, one of these length scales becomes of 
order L, the long time dynamics of the system is totally controlled by the artificial 
periodic boundary conditions, and the final states are artifacts of these unphysical 
boundaries. 

The above discussions show that the homogeneous and inhomogeneous shearing 
states, as well^alhe kinetic state, discussed in the previous sections, and reported in 
the literature,Hliy'til do not represent physical properties of macroscopic (idealized) 
granular flows, but are caused by the artificial periodic boundary conditions. It 
is also interesting to observe that in cases where the boundary conditions have 
been modelled in a physically more realistic way and in different geometriest£l — 
i.e. by a wall kept at a constant temperature — different asymptotic states have 
been observed. This discussion suggests that more realistic modelling of collisions 
between granular particles and physical walls is imperative for understanding the 
asymptotic states of granular flows contained in finite geometries. 
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